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ABSTRACT 



Thermal conduction has been suggested as a possible mechanism by which suffi- 
cient energy is supplied to the central regions of galaxy clusters to balance the effect 
of radiative cooling. Recent high resolution observations of the nearby Virgo cluster 
make it an ideal subject for an attempt to reproduce the properties of the cluster by 
numerical simulations, since most of the defining parameters are comparatively well 
^j- \ known. Here we present the results of a simulated high-resolution, 3-d Virgo cluster 

■ for different values of thermal conductivity (1, 1/10, 1/100, times the full Spitzer 

c*^) value) . Starting from an initially isothermal cluster atmosphere we allow the cluster to 

I/"") , evolve freely over timescales of roughly 1.3 — 4.7 x 10 9 yrs. Our results show that ther- 

' nial conductivity at the Spitzer value can increase the central ICM radiative cooling 

time by a factor of roughly 3.6. In addition, for larger values of thermal conductvity 
the simulated temperature and density profiles match the observations significantly 
I 1 better than for the lower values. However, no physically meaningful value of thermal 

conductivity was able to postpone the cooling catastrophe (characterised by a rapid 
increase in the central density) for longer than a fraction of the Hubble time nor ex- 
plain the absence of a strong cooling flow in the Virgo cluster today. We also calculate 
the effective adiabatic index of the cluster gas for both simulation and observational 
data and compare the values with theoretical expectations. Using this method it ap- 
pears that the Virgo cluster is being heated in the cluster centre by a mechanism other 
than thermal conductivity. Based on this and our simulations it is also likely that the 
thermal conductvity is suppressed by a factor of at least 10 and probably more. Thus, 
we suggest that thermal conductvity, if present at all, has the effect of slowing down 
the evolution of the ICM, by radiative cooling, but only by a factor of a few. 

Key words: thermal conduction-galaxies: clusters-individual(Virgo)-galaxies: active- 
cooling flows 



1 INTRODUCTION 

The strong central peak of the X-ray surface brightness 
profiles of many clusters of galaxies is generally inter- 
preted as the signa ture of a cooling flow llCowie fc McKeel 
Il977l: lFabiar] ll994h Standard cooling flow theory suggests 
that this peak should coincide with the deposition of large 
amounts of cold gas and significant excess of star forma- 
tion at the centre of the cooling flow. However, the XMM- 
Newton and Chandra satellites have demonstrated both the 
lack of the e xpected cold gas, at temperatures below roughly 
1 keV (e.g. lEded 1200 J) , and that spectroscopically mea- 
sured mass-deposition rates are roughly one tenth of the 
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value inferr ed from the peaks of clu ster X-ray surface bright- 
ness maps JVoiet fc Fabianll2004) . These findings suggest 
that the gas in the cluster cores is reheated; either con- 
tinually or periodically, in order to produce the observed 
minimum central temperature (or entropy) and the lack 
of significant excess star formation. In recent years, work 
has concentrated on two main possibilities for heating the 
gas in the central r egions of clusters: (1) heating by out- 
flows from AGN (e.g.lTabor fc Binnevll993tlChurazov et alJ 
l200ll iBriiggen fc Kaiserl 120021)" and (2) transportation of 
heat from th e outer regions of the cluster by thermal con- 
duction (e.g. IZakamska fc Naravanll2003l : IVoigt et al1l2002l : 
IVoigt fc F abian 200 J). In this paper we shall concentrate 
purely on thermal conduction. 

The presence of large temperature gradients in the 
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centres of many relaxed clusters allows the possibility 
that thermal conduction may play a role in reducing 
the radiative energy l osses in the cluster centre. Indeed, 
IZakamska fc Naravanl i2003T) show that for several galaxy 
clusters the inward heatflow from the outer regions due to 
thermal conduction may be sufficient to balance the radia- 
tive losses although, in a different sample IVoigt et alJ <2002h 
show that thermal conduction is only able to balance radia- 
tive losses in the outer regions of the cluster core. However, 
though it is possible to achieve an energy balance for these 
clusters in their present state, this requires the thermal con- 
ductivity to vary as a function of radius. Furthermore, it is 
unclear what effect thermal conduction would have on the 
intra cluster medium (ICM) if it constitutes a significant 
process throughout the lifetime of the cluster. Since thermal 
conduction acts to oppose the formation of the temperature 
gradient that causes it, one might expect that the temper- 
ature and density profiles for those clusters in which ther- 
mal conduction plays a significant role to appear different 
compared to those clusters in which thermal conduction is 
less effective. In particular, the temperature profiles in such 
clusters could be much flatter than in clusters where thermal 
condcution is heavily suppressed. In addition, temperature 
gradients also occur outwards towards the edge of clusters. 
Therefore, it is possible that the high thermal conductivi- 
ties required to balance the radiative losses may also have a 
cooling effect by transporting energy away from t he cluster 
centre and out into intercluster space dLoebll2002l) . 

One of the primary problems for studying the effects of 
thermal conduction in plasmas is the unknown value of the 
thermal conductivity. The theoretical value for the thermal 
conductivity of a fully ionized, unmag netized plasma was 
originally calculated by ISpitzerl fl962). However, observa- 
tions indicate t he existence of a magn etic field in the ICM 
(B < 10 fj,G; e.g. lCarilli fc Tavlorl2002T) which can affect this 
value dramatically. For magnetic fields of these strengths the 
electron and proton gyro radii are much smaller than their 
mean free paths along the field lines. The precise effect that 
magnetic fields have on thermal conductivity is unclear, al- 
though for magnetic fields of the observed strengths the ef- 
fective ther mal conductivity is dete rmined by the topology 
of the field iMarkevitch et al"ll2003l) . A magnetic field tan- 
gled on scales much longer than the electron mean free path 
will suppress thermal conduction perpendicular to the field 
lines reducing the thermal conductivity to at most one third 
of the full Spitzer value ( if transport along the field line is 
unaltered) iSarazin|ll98rf) . For a magnetic field tangled on 
scales shorter than the electron mean free path the conduc- 
tivity is suppressed by a factor of up to 100 as the electrons 
have to travel furt her to diffuse a given distance perpen- 
dicular to the field foibbldll989T) . In addition, for tangling 
lengths that are comparable to the electron mean free path, 
conduction along the magnetic field line can be reduced by 
a factor of roughly 10 by the effect of magnetic mirrors 
iMalvshkin fc K ulsrudll2 00lF). Alternatively , rece nt theoret- 
ical w ork by iNaravanfe Medvedevl teOOll) and ICho et all 
( 2003) has shown that for certain spectra of field fluctu- 
ations, turbulent magnetic fields are less efficient at sup- 
pressing thermal conduction than previously thought. The 
highly tangled, turbulent magnetic field may allow signifi- 
cant cross field diffusion such that the thermal conductivity 
remains of the order of the Spitzer value. In any case the 



collisional thermal conductivity on 100 kpc scales is likel y 
to be reduced by a factor of 3-100 iMarkevitch et aTl feoOS"). 

We define the suppression factor, /, by K = /ts, where 
k is the actual thermal conductivity and «s is the Spitzer 
thermal conductivity. To date this fa ctor has been very diffi- 
cult to constrain from observations: IVoigt et alJ (120021) and 
IZakamska fc Naravanl J2003T) find that suppression factors of 
> 0.3 are sufficient to balance radiative losses in some cluster 
while other cluster require unphysically high values (/ > 1). 
In othe r cases, the temperatu re drops across cold fronts in 
A2142 fettori fc Fabianl 20001 implies a suppression factor 
of < 0.004. At the boundaries be tween the ICM and ga laxy- 
size dense gas clouds in Coma dVikhlinin et al 1 l200ll) the 
thermal conductivity is found to be suppressed by a factor 
of order one hundred. However, one might reasonably expect 
both of these cases not to be characteristic of the bulk ICM, 
since they contain boundaries between different gas phases 
and disjointed magnetic fields. Another est imate uses the ex- 
isten ce of cold filaments in galaxy clusters iNipoti fc Binnevl 
2004) to put constraints on the value of the suppression fac- 
tor. Although not all of the values of the required quantities 
are well known, they find that for the Perseus cluster the 
suppresion factor is < 0.04, whic h result they cons i der to 
be compatible with the results of Ma rkevitch et all J2003T) 
based on observations of temperature gradients present in 
A754. 

The aim of the work presented here is to investigate 
whether thermal conduction can prevent the catastrophic 
radiative cooling of the gas at the centres of galaxy clusters 
by transporting thermal energy from the cluster outskirts 
to the centre. We also investigate whether this process gives 
rise to gas temperature and density profiles compatible with 
those derived from X-ray observations. 

The plan of this paper is as follows. In Section|5|we dis- 
cuss why thermal conduction would be an unlikely solution 
to the cooling flow problem from analytical predictions. Sec- 
tion describes the numerical model and code which we 
use to obtain our results. Section ^. 2l states the initial condi- 
tions specific to the problem we are investigating. Section HPl 
gives the details of all simulations, e.g. time simulated and 
value of thermal conductivity used. In Section 4 we present 
the results of our simulations in the form of temperature and 
density profiles, emissivity profiles, the effective adiabatic in- 
dex, central temperature and density, and mass flow rates 
respectively. We also discuss the comparisons of these re- 
sults with recent observations of the Virgo cluster and their 
implications. We present our conclusions in Section |S] and 
demonstrate the numerical convergence of our results in the 
Appendix. 



2 HEATING BY THERMAL CONDUCTION 

Thermal conduction transfers heat so as to oppose the tem- 
perature gradient which causes the transfer. Therefore, if 
thermal conduction occurs at all in galaxy clusters, it will 
certainly result in some heating of the cluster centre, but 
to what extent? If thermal conduction is to provide a solu- 
tion to the cooling flow problem it must fulfill the following 
conditions: 

(i) temperature and density profiles must be simultane- 
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ously comparable with observation s and stable for cosmo- 
logically significant timescales (e.g. lAllen et al. 2001). 

(ii) the deposition of large amount of co ld gas in th e cen- 
tre of the cluster must be suppressed (e.g. Edge 2001). 

In order to assess the likelyhood that thermal conduc- 
tion is capable of satisfying the above requirements we make 
analytical predictions based on the relevant hydrodynamic 
equations. 

The rate of radiative energy loss per unit volume 
(e ra d) is given by equation Q where n c is the electron 
concentration, n p is the concentration of protons/positive 
ions and A ra d (T, Z) is the cooling function 



where A rad = 2.1 x 10~ 27 VT erg s^cm" 3 
iRvbicki fc Lightmanl Il979l) is the cooling function for 
pure bremsstrahlung emission, T is the gas temperature 
and n e is the electron number density. 

Using equation J7J, with the initial conditions near the 
centre of the Virgo cluster (see section l3~2t . the central cool- 
ing time is approximately 200 Myr. This is comparable to 
the results of our simulation without thermal conduction. 
The cooling time at ~ 100 kpc is ~ 10 12 yrs. 

Given the form of equation , the central cooling time 
for a radiatively cooling and thermally conducting ICM is 
given by 



trad = ~n e n p A ra d (T, Z) [erg s 'cm 



(1) 



The rate of heating due to thermal conduction is given 
by the divergence of the thermal flux 
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where k is the thermal conductivity defined in the introduc- 
tion. 

For the high temperatures of cluster atmospheres the 
coefficients of thermal co nductivity and viscosity are given 
by the value calulated by ISpitzerl jl962T) . In the absence of 
magnetic fields the coefficients of Spitzer thermal conduc- 
tivity and viscosity are given by 
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(4) 



where lnA c is the Coulomb logarithm and A c is given by 
equation 0. 

Due to the large variation in density over cluster-scale 
distances it is necessary to take into ac count the variatio n 
of the Coulomb logarithm given by re.g. lChoudhHrilll99^ 



^cool ~ [SJ 
^rad ^cond 



(8) 



where e ra d and e con d are defined by equations Q and (|5J 
respectively. 

Thus, if thermal conductivity is the sole process by 
which the radiative cooling is opposed, then the predicted 
cooling time of 200 Myr, based on observations, must be ex- 
tended by at least an order of magnitude. From equation (|HJ 
this requires the energy flux due to thermal conduction to 
be equal to the radiative energy flux to within better than 1 
%. Such a situation may be possible for exceptional clusters 
e.g. those with very low density, so that the radiative losses 
are small, and high temperatures, so that the thermal con- 
ductivity is large. However, due to the independent nature 
of the competing processes it seems unfeasible to expect this 
to be the case in every galaxy cluster. In addition, although 
thermal conductivity probably increases the cooling time of 
the gas at the cluster centre by a factor of a few, in doing so 
it also reduces the cooling time at larger radii introducing a 
secondary cooling problem. 

However, we note that the above analysis is based only 
on estimates, since it is impossible to solve the full hydrody- 
namic equations analytically. Therefore, in order to better 
understand the effect of thermal conduction on the ICM we 
perform 3-d simulations in which the hydrodynamic equa- 
tions are solved numerically. 



A c = 247T71, 



/ 87re 2 n c 
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The opposing energy fluxes of radiative cooling and 
thermal conduction determine the behaviour of the cluster 
gas. For a spherically symmetric cooling flow, in which we 
ignore the effect of viscous heating and gravitational effects, 
the energy equation is 



d / 3fcBfleT 

dt \ 2/im p 



-Crad+Econd [ergS 1 Cm 3 ] 



(6) 



where the symbols have the definitions given above. 

A rough estimate of the cooling lifetimes of the ICM 
can be used to illustrate the point that thermal conduction 
would be an unlikely solution to the cooling flow problem. 
For the case in which thermal conduction is completely su- 
pressed the cooling time of the gas can be defined as 



(7) 



3 NUMERICAL MODEL 



3.1 The Code 



We performed high resolution, 3-d numerical hydrodynamic 
simulations to produce temperature and density profiles of 
model Virgo clusters. For simplicity we assume that the dark 
matter halo of our simulated cluster is fully formed and 
stationary throughout the simulation. We also neglect the 
self-gravity of the gas. Thus the gravitational potential con- 
fining the gas is fixed. We also self-consistently include the 
effects of radiative cooling, thermal conductivity and vis- 
cosity. We do not include any external heating mechanisms 
such as AGN. By neglecting the self-gravity of the cluster 
gas we are able to run the simulations over cosmologically 
relevant timescales (few x 10 8 — few x 10 9 yrs). 

Our simulations were performed with FLASH2.3, an 
Adaptive Mesh Refinement (AMR) hydrodynamical code 
developed and ma de public by the AS CII Center at the Uni- 
versity of Chicago jFrvxell et al.l2 0001. FLASH is a modular 
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block-structured AMR code, parallelised using the Message 
Passing Interface (MPI) library. 

FLASH solves the Riemann problem on a carte- 
sian grid using the Piec ewise-Parabolic Method (PPM; 
I Woodward fc Colellalll984h . It uses a criterion based on the 
dimensionless second derivative D 2 = Fd 2 F/dx 2 / (dF/dx) 2 
of a fluid variable F to increase the resolution adaptively 
whenever D 2 > ci and de-refine the grid when D 2 < a, 
Ci,2 are tolerance parameters. When a region requires refin- 
ing, child cells of half the size of the parent cells are placed 
over the offending region, and the coarse solution is inter- 
polated. In our simulations we have used density, pressure 
and temperature as the re finement fluid variable F (see also 
iDalla Vecchia et aljEooll . 

FLASH imposes the analytic gravitational potential to 
the grid cell, and computes the corresponding gravitational 
acceleration. This is done only once since we neglect the 
self-gravity of the gas. To interpolate the intial density field 
to the FLASH mesh, we impose a uniform initial grid with 
a sphere of higher resolution in the central regions with a 
radius of 16.2 kpc, because of the steeper gradients in the 
cluster centre. This is the minimum refinement allowed dur- 
ing the simulation, the maximum refinement was set to 9 re- 
finement levels, equivalent to 2048 3 uniform computational 
cells. We use periodic boundary conditions and a computa- 
tional box size large enough (648 kpc) that the density at 
the cluster edge remains approximately constant over the 
lifetime of the simulation. 

FLASH calculates three limiting timesteps based on the 
requirements for consistency in the hydrodynamics, radia- 
tive cooling and diffusive processes, i.e. thermal conduction 
and viscosity, of the solution in the entire computational 
grid. The solution is then advanced by choosing the mini- 
mum of these timesteps. The hydro timestep is derived by 
the Courant condition dt — CSx/cs where Sx is the cell di- 
mension, c s is the sound speed in that cell and C is the CFL 
coeffient which is set to 0.8 in our simulations. The cool- 
ing timestep is given by: dt = ei/(dei/di) where ei is the 
internal energy per unit volume and dei/dt is the rate of 
change of internal energy. We usually find this timestep to 
be very large except for very high cooling rates. The diffu- 
sion timestep is given by: dt = 0.5(min(da; 2 ))/max(x) where 
min(da; 2 ) is the smallest cell dimension, x ls the diffusivity 
and max(x) is the largest diffusivity in the computational 
zone. This timestep is smaller for larger values of thermal 
conduction and viscosity. 

The simulations were performed using Iridis, one of the 
University of Southampton's Beowulf clusters. Iridis has 115 
computational nodes connected with myrinet of which we 
usually used 8 to 10 nodes (16 to 20 processors). 

3.2 Initial Conditions 

From the assumption that the Virgo cluster is currently ap- 
proximately in hydrostatic equilibrium it is possible to de- 
rive the gravitational acceleration as a function of radius 
within the cluster if both the temperature and the density 
distributions of the gas are well defined. For these we use 
the best-fit functions for the te mperature and el e ctron num- 
ber density distributions from IGhizzardi et alj i2004l) who 
combine the data from Chandra, XMM-Newton and Beppo- 
SAX. They fit the deprojected temperature and electron 



number density profiles with a Gaussian (see equation [5J 
and double beta profile (see equation I1UI respectively. We 
consider this method more accurate than estimating the ap- 
propriate parameters for a Navarro, Frenk & White potential 
iNavarro et alll 997l. since the gravitational potential in the 
central regions of a real galaxy cluster will be dominated by 
the central galaxy. 



7 = /;-7M-q>( -~ 



K], 



the best fit values are T = 2.78 xlO 7 K , Ti -■ 
K and a = 7.39 xlO 22 cm. 

The electron density is given by: 

_ no ni r 

Uc ~ (1 + (r/ro) 2 ) a ° (1 + (r/n) 2 )"i [ ' 

the best fit values are given by no = 0.089 cm" 



(9) 



8.997 x 10 f) 



(10) 



1.6 xl0 22 cm, 



7.39 xlO 22 cm 



m = 0.019 
a = 1.52, 



cm , r 
Qi = 0.705. 

The double /3-pr ofile for the gas density found by 
IGhizzardi et al.l (120041) is interesting given the observation 
that although in general cluster temperature profiles are con- 
sistent with a single phase gas, the X-ray surface brightness 
is less centrally peaked than expected. This is taken as evi- 
dence for distributed mass d eposition througho ut the flow, 
caused by a multiphase gas jFabian et al.l f2002. and refer- 
ences therein). 

Work by lEvrardl Jl99fJ) and iNavarro etaD dl995l) 
showed that the central megaparsec of clusters are roughly 
isothermal due to the effect of shock heating of the infall of 
the gas in the formation of the cluster. Therefore our initial 
conditions assume that the ICM is uniformly heated to the 
cluster Virial temperature and has an initial density profile 
that ensures hydrostatic equilibrium with the gravitational 
potential derived above. The assumption of a precisely uni- 
form initial temperature is adequate for the purposes of this 
work since we simulate only the central few hundred kilopar- 
secs of the Virgo cluster. However, we note that it may still 
be simplisitic with regard to residual temperature and den- 
sity perturbations arising from the formation of the cluster 
which should be taken in to account for a complete study 
of ICM evolution. However, this will be dealt with in future 
work. 

We estimate the Virial temperature (T v i r ) of the Virgo 
cluster from 



GMvirfirrip 

kBRvii 



[K], 



(11) 



where M v i r is the Virial mass, 7? v i r is the Virial radius, G is 
the Universal gravitational constant, and fj,m p is the mean 
molecular mass of the gas. 

We are able to estimate the Virial mas s and radius us- 
ing th e gravitational mass profile given bv IGhizzardi et 

ID 

(2004). We adopt a Virial mass of 10 13 M Q with a Virial ra- 
dius of 100 kpc and find the Virial temperature of the Virgo 
cluster to be ~ 3.1 x 10 7 K. The initial value for the uniform 
temperature was determined by finding a temperature near 
to the Virial temperature, for which the simulated tempera- 
ture profile, after a period of 10 s yrs of evolution of the ICM, 
was comparable to current observations. We found that an 
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Figure 1. Cooling function as a function of plasma temperature. 
Note that the effects of line cooling at lower temperatures mean 
that the gas cools much more rapidly at lower temperatures than 
if we just assume bremsstrahlung. At high temperatures the cool- 
ing is dominated by bremsstrahlung. 



initial temperature of 3 x 10 7 K resulted in an adequate 
temperature profile. 

The initial density profile required for hydrostatic equi- 
librium is determined by the derived gravitational potential 
and the initially uniform temperature up to a multiplica- 
tive constant. This constant is defined as the density at the 
cluster centre. Since the gas density in the cluster outskirts 
will not vary much throughout the simulated time, we ad- 
just this constant such that the density in the outer regions 
of our cluster is comparable to the currently observed den- 
sity in these localities. Thus we found an appropriate central 
density of 5.15 x 10~ 26 g cm -3 . 

In our model, the ICM is a single, ideal fluid, of half- 
solar metallicity which is assumed to be fully ionised with 
an adiabatic index of 5/3. For such a composition the total 
mass density of the ICM is given by: p = (8/7)m p n c where 
n e is the electron number density and m p is the proton mass. 

In order to properly take into account the effect 
of chemical enrichment and metallicity on the radia- 
tive cooling we adop t the cooling function calculated by 
ISutherland fc Dopital <ll993t) . The function is suitable for 
low density astrophysical plasma and was obtained from 
the study of plasmas cooling in both equilibrium and non- 
equilibrium situations in the temperature range 10 4 - IO 85 
K. By using this cooling function we are able to account 
for the effects of emission-line cooling which become signifi- 
cant for temperatures below ~ 10 7 (see figure^- We assume 
the metallicity to be half of the solar value fedee fc Steward 
Il99ll) and that the cooling function depends only on tem- 
perature and metallicity, but not density. For temperatures 
below 10 4 K the cooling effect is 'switched off' meaning that 
we do not accurately model cooling below this temperature. 
However, for none of our simulations were such low gas tem- 
peratures reached. 

In order to implement the plasma thermal conductivity 
correctly it is essential to know whether the electron mean 
free path is less than the scale length of the temperature gra- 
dient. The scale lengt h of the temp erature gradient can be 
defined as h = T/VT {g arazinll98€f) . For electron mean free 



paths which are greater than the scale length of the tempera- 
ture gradient the thermal conduction is said to 'saturate' and 
the h eat flux approaches a limiting value llCowie fc McKeel 
Il977l) . For electron mean free paths much less than the tem- 
perature gradient the heat flux depends on the coefficient of 
thermal conductivity and the temperature gradient. 

The mean free path of an electron (A e ) in an unmagne- 
tized plasma is given by dSarazinlll986^ 

A c = - 1/2 V [cm] (12) 

47T 1 /2 noe 4 m A c 

Substituting typical values for near the centre of the 
cluster (n c ~ 0.1cm~ 3 ,T c ~ 10 7 K) gives an electron mean 
free path of roughly 10 pc where the length of the temper- 
ature gradient is of the order of kiloparsecs. For the outer 
cluster regions (n e ~ l(T 4 cm- 3 ,T c ~3x 10 7 K) we find 
a mean free path of ~ 20 kpc where the scale length of 
the temperature gradient is of the order of hundreds of kpc. 
This means that we are able to use the usual, 'non-saturated' 
form for the thermal conductivity at all radii for the Virgo 
cluster. 

Using the above result, we implemented thermal con- 
duction and viscosity with the transport coefficients given 
in equations (3| and gj, suppressed by the factor, /, into 
the simulation software. We use the same suppression fac- 
tor for both conductivity a nd viscosity as the two processes 
are intimately linked (e.g. iKaiser et alJl2005T) . In addition 
we also include a control run in which thermal conduction 
and viscosity are absent. For the cases in which we consider 
the effects of suppressed thermal conductivity and viscosity, 
the suppression factor is constant throughout the ICM. In 
reality the suppression factor is unlikely to be uniform due 
to varying magnetic field strength and tangling length. We 
note that unless the central regions of galaxy clusters are 
turbulent such that the t hermal conductivity is supp ressed 
only by a factor of a few iZakam ska fc Naravanll200Stl , then 
it is likely that the tangling length of the magnetic field is 
shorter in the central regions (due to compression). There- 
fore the thermal conductivity may be smaller in the centre 
compared to the outer regions. The effect of this would be 
to reduce thermal conduction where it is required most. 

When taking into account thermal conductivity as a 
source of heat for balancing radiative losses, i.e. e con d ~ trad, 
it is necessary to have sufficient energy in the outer regions 
of the cluster to achieve this. We can consider the outer 
regions of a cluster to be an infinite heat bath, however, 
larger volumes take longer to simulate. Simulating smaller 
volumes is quicker, but it may starve the simulated cluster 
of the energy that thermal conduction requires. The size of 
the computational box was fixed as a cube of dimension 648 
kpc, with the cluster centre situated at the box centre. For a 
computational box of this size the regions outside 5 kpc (in- 
side which the cooling catastrophe was observed to occur in 
preliminary simulations) contain roughly 10,000 times more 
energy than can be supplied to the central 5 kpc over the 
lifetime of a simulation (roughly 2-5 Gyr). 

It is also wort h pointing out that a problem noted by 
iDolae et alJ i2004T) in which simulations involving radiatve 
cooling and star formation found an increase in tempera- 
ture at the cluster centre was also observed in our prelim- 
inary simulations. We attribute this effect to limited nu- 
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merical resolution, since this effect of compressional heating 
was eliminated by increasing the number of computational 
blocks in the central region of our cluster (see Appendix). 



3.3 The Simulations 

The nature of the problem we are considering here is 1- 
dimensional in many respects; this includes the initial con- 
ditions and the non-evolving gravitational potential. In gen- 
eral the physical mechanisms operating in this setup will re- 
sult in mainly radial flow of the gas. It is possible to achieve 
much higher spatial resolution in a 1-d compared to a 3-d 
simulation. However, by employing a 3-d geometry we allow 
for any non-radial processes which may occur, for example: 
the growth of non-spherical modes of the thermal instabil- 
ity or non-radial mass flow which is likely to occur near the 
cluster centre. 

In order to determine any differences between employ- 
ing 3-d or 1-d coordinate systems we have performed 1-d 
spherically symmetric simulations for two cases: i) the ab- 
sence of thermal conduction and viscosity and ii) the limit of 
full Spitzer thermal conductivity and viscosity. The results 
of these simulations are similar to those performed using 
a 3-d coordinate system, and are given more completely in 
Appendix^] The main observable difference is that the evo- 
lution of 1-d spherically symmetric clusters is slightly more 
rapid than in the 3-d case. We consider this behaviour con- 
sistent with the inherent differences between the hydrody- 
namic equations in 1-d and 3-d. For example, there are extra 
dissipation terms both including and excluding the viscos- 
ity in the 3-d case which are absent in the 1-d case (see 
Appendix [5J thus any non-spherically symmetric processes 
will never be allowed to develop in 1-d. We therefore proceed 
with 3-d simulations since we believe that the 3-d effects are 
important in the overall evolution of the gas. 

We performed 4 simulations to which we assign labels 
such that the number preceeding n is the suppression factor, 
/• 

(i) radiative cooling alone (zero thermal conduction and 
viscosity) (Ok). 

(ii) radiative and thermal conductivity at one hundredth 
of the full Spitzer value (0.01k). 

(iii) radiative and thermal conductivity at one tenth of 
the full Spitzer value (0.1k). 

(iv) radiative and thermal conductivity at the full Spitzer 
value (Ik). 

Table 1 provides a summary of the simulations. 

Each simulation was stopped after the onset of a cooling 
catastrophe in the cluster centre. As a criterion for a cooling 
catastrophe we used the timestep adopted by the simulation 
software. If this timestep became insignificant compared to 
the duration of the entire simulation, we stopped the sim- 
ulation. This occured for central temperatures of less than 
~ 10 7 K in all simulations. 



4 RESULTS 

4.1 Temperature and Electron Number Density 
Profiles 

We present spherically averaged profiles, from our simula- 
tions, for both temperature and density, rather than 1-d 
slices through the cluster. This is done by defining a num- 
ber of spherical concentric shells at a number of radii, in our 
case 200, from the cluster centre and averaging the tempera- 
ture and density between adjacent shells. The temperatures 
and densities are compared to observations in figures [5] to [5] 

The temperature and electron number density profiles 
are plotted at intervals of 3.17 x 10 s yrs for all, except simu- 
lation Ik which is plotted for every 6.34 x 10 s yrs, and for the 
final time of the simulation (see table 1). We also include the 
commen initial conditions for simulation 0k only. For each 
simulation we note that the central density increases at all 
times during the ICM evolution. However, the behaviour 
of the gas temperature is rather more complex: during the 
early phases of evolution of the ICM the temperature falls, 
but appears to rise after the onset of the cooling catastrophe 
(see also section HO} . 

The cooling catastrophe is characterised by the rela- 
tively sudden and dramatic increase in the central density 
over a short period of time (few x 10 8 yrs) . It should be 
noted that no physical value of thermal conductivity can 
prevent the occurance of a cooling catastrophe. 

One result from our simulations is that regardless of 
whether thermal conductivity and viscosity are present, 
both the temperature and density profiles eventually ap- 
proach generic profiles. After 2 x 10 9 yrs even the simula- 
tion with full Spitzer thermal conductivity develops into a 
cooling catastrophe with temperature and density profiles 
similar to the other simulations. 

The temperature profiles are characterised by a narrow, 
but deep central dip which deepens with time such that the 
central temperatures fall below the observ ed approximate 
minimum of T vir /3 (e.g. lAllen et al.ll200ll) . Qualitatively, 
the main effect of thermal conduction is an increase of the 
width of the central temperature dip. This is due to ther- 
mal conduction tranporting energy down the temperature 
gradient towards the cluster centre. This reduces the rate at 
which the temperature falls in the cluster centre, but also 
increases the effective rate at which the gas cools at larger 
radii. The time taken to evolve into this profile is roughly 
inversely proportional to the magnitude of the thermal con- 
ductivity. 

The broader temperature dip associated with larger 
thermal conductivities is accompanied by larger densities in 
the same region of the cluster. This has the same physical 
explanation as above. Due to thermal conductivity, the gas 
at larger radii loses more energy than it would otherwise do 
due to radiative losses. This prompts the inflow of gas at 
larger radii creating a build-up of material in a larger region 
of the cluster centre than for lower values of the thermal 
conductivity. 

The density profiles can be adequately described by 
a strongly centrally peaked distribution in which the peak 
grows with time. As energy is radiated away most rapidly in 
the densest regions, the pressure support is lost from these 
central regions and the weight of the overlying layers com- 
presses the central gas thus increasing the density dramati- 
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Name 


Suprcssion factor 


Simulated time (yrs) 


Ok 





1.29 xlO 9 


O.Olft 


1/100 


1.35 xlO 9 


0.1k 


1/10 


1.89 xlO 9 


Ik 


1 


4.72 xlO 9 



Table 1. Summary of the five main simulations 




distance [kpc] distance [kpc] 



Figure 2. Temperature and density profiles evolving with time for simulation 0k (see table 1). The thick lines for both temperature and 
density arc the functions fitted to the data points (triangles) by Ghizzardi et al. (2004). The top line in the temperature plot shows the 
temperature profile after 3.17 X 10 8 yr and the bottom line at time of the end of the simulation. These times are given in table 1 for 
each simulation. The intermediate lines represent the temperatures at intervals of 3.17 X 10 8 yr after the top temperature profile. The 
temporal sequence of the lines is reversed (bottom to top) in the density plot. 
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Figure 3. Temperature and density profiles evolving with time for simulation 0.01k. 
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Figure 4. Temperature and density profiles evolving with time for simulation 0.1k. 




Figure 5. Temperature and density profiles evolving with time for simulation Ik. Note that in this figure the temperature and density 
are plotted at intervals of 6.34 X 10 s yrs. 



cally. The fall in central temperature and rise in central den- 
sity that we observe are both the obvious signs of a strong 
cooling flow and eventually a cooling catastrophe in which 
large quantities of gas are being deposited into the central 
regions of the cluster. 

It is evident from figures|2]to|S]that for increasing ther- 
mal conductivity the fit of the simulated density profiles to 
the observational data improves at larger radii. However, 
the profiles are not stable over sufficiently long timescales, 
e.g. the Hubble time. Thus, thermal conduction only delays 
the occurance of a cooling catastrophe and cannot simula- 
taneously reproduce temperature and density profiles that 
are consistent with observations. 

4.2 Emissivity Profiles 

From the simulated density and temperature data it is possi- 
ble to derive X-ray surface brightness maps for our radiative 



cooling function. Previously we have compared the density 
and temperature distributions of our simulated clusters with 
those derived from observations of the X-ray surface bright- 
ness of the gas, since these are the fundmental variables used 
by FLASH. However, calculating the surface brightness for 
our simulations allows for a more direct comparison with the 
observational data. In practice the surface brightness along 
a given line of sight is dominated by emission from gas at the 
smallest radii from the cluster centre along that particular 
line of sight. Thus we calculate the emissivity of the gas (e), 
i.e. the rate of energy loss per unit volume, and plot this as 
a function of radius. The bolometric emissivity is given by 

e = n e n p A rad (T, Z) [erg s^cm" 3 ] (13) 

Figures |S| to [7] show the emissivity after 6.34 x 10 8 yrs 
and at the final point in the simulation compared to the 
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currently observed em issivity calculated from the work of 
iGhizzardi et~ail fo004\ . 

XMM-Newton observations revealed that the X-ray sur- 
face brightness wa s less centrally peaked than expected (e.g. 
Fabi an et al . 2002). From figurcsUJtoQit is evident that our 
simulated clusters always develop a strong central peak in 
the surface brightness distribution, whether or not thermal 
conduction is taken into account. Such a peak is to be ex- 
pected given the density profiles in section 14.11 but is not 
present in the observational data. This again indicates that 
an extra heating mechanism in the central few kiloparsecs 
is required to prevent the excessive central cooling. 

We note that the general fit at large radii to the obser- 
vational data seems to improve as the value of the thermal 
conductivity is increased. In addition, the emissivity pro- 
files show more excess emission out to larger radii for the 
simulations with a larger thermal conductivity. Thus, the 
emissivity plots re-iterate the previous observation from the 
density profiles that for higher thermal conductivities there 
is more mass in the central 20-50 kpc than for lower values 
of thermal conductivity. This is consistent with the mass 
flow rates which are larger at all radii for later times in 
the presence of high thermal conductivity than those for 
lower thermal conductivity (see section l4.5[l . Hence, in or- 
der to better match the observations by preventing a cooling 
catastrophe it appears that the presence of a larger thermal 
conductivity actually requires the input of more thermal en- 
ergy. In addition, this heat source would need to distribute 
its energy out to larger radii for the cases of larger thermal 
conductvity. The extra energy input required is roughly pro- 
portional to the extra volume which requires heating, and 
so an ICM with large thermal conductivity could require 
100-1000 times more energy than for smaller thermal con- 
ductivities. 



4.3 Effective Adiabatic Index 

As an additional diagnostic for the behaviour of the clus- 
ter gas we calculate the effective adiabatic index defined by 
7 c ff = dlnP/dlnp where P is the gas pressure, and p is 
the ga s density. Following the method of lBirnboim fc DekeJ 
(2003) we calculate 7 e fi for a situation in which the volume 
of gas is constant, but is subjected to heating and cooling. 
The pressure in terms of the internal energy per unit mass, 
e, is given by 

P = (7 - l)ep [erg cm -3 ], (14) 

where 7 is the true adiabatic index (e.g. 5/3 for a monatomic 
gas) 

The rate of change of gas pressure with time is given by 

P = (7 — l)(ep + ep) [erg cm -3 s~ ]. (15) 

The rate of change of internal energy per unit mass for 
a constant volume of gas subjected to heating and cooling 
is given by 

e = h-q [erg g _1 ], (16) 



where h and q are the heating and cooling functions per unit 
mass respectively. 

By converting all the relevant quantities to be compat- 
ible with those we have defined previously and using the 
definition of the effective adiabatic index we find 7 e g to be 
given by 

Taff = l + ( 7 -l)^ ^ j, (17) 

where H is the heating function (the heating analogue of the 
cooling function, A ra d), h is the rate of change of number 
density with time and the other symbols have their usual 
definitions. 

In this derivation the nature of the heating and cool- 
ing mechanisms is not important, however, we choose to 
represent them in terms of quantities that we have already 
defined. From equation 11711 we see that in the limit where 
either there is no heating or cooling or these competing ef- 
fects balance, the effective index is 1, which corresponds to 
an isothermal equation of state. Furthermore, we note that 
in the cases where no heating occurs, j e g is always less than 
one and for extreme cases will tend to zero as h tends to 
ra/icooi (where t coo i is the radiative cooling time). For a con- 
stant pressure cooling flow f e s will be equal to 1/3. In the 
cases where heating occurs this will have the effect of increas- 
ing 7eff. Using this equation we predict for the simulations 
which include thermal conductivity at a significant fraction 
of the Spitzer value that the gas exhibits a 7 c fj that is greater 
than for lower values of thermal conductivity, at least in the 
cluster centre where radiative cooling is strongest. 

In figure [S] we compar e the effective ad iabatic index ob- 
tained from the results of IGhizzardi et al] (120041) with 7 e g 
for simulations 0k and Ik. We note that in the absence of 
thermal conductivity j e g is roughly equal to 1 (isothermal) 
in our simulation only in the outer regions of the cluster, 
but falls steeply inside the central 10 kpc towards the cen- 
tral value given in table 2. The upturn at roughly 3 kpc is 
probably due to compressional heating towards the end of 
the simulated time. We also note that there is a minimum 
at roughly 25 kpc and a peak at roughly 15 kpc. For the full 
Spitzer value of thermal conductivity, y e g follows roughly 
a similar pattern to the above except that the variation is 
smoother. In addition, the central effective adiabatic index 
is higher due to the heating effect of thermal conductivity. In 
comparison, it is evid ent that for the functio ns fitted to the 
observational data bv IGhizzardi et al] ||2004|) 7 e fj is isother- 
mal throughout the cluster except for around 30 kpc where 
we observe a minimum similar to that noted in the simula- 
tion data. Of particular interest is the fact that there are 4 
data points between the radii of 4-12 kpc with an effective 
adiabatic index of roughly 1 that are not consistent with 
either of the two simulations. By equation I17H this value of 
the observed y e g must be due to heating balancing cooling. 
From the simulations we can see that such a profile for j c s 
is unattainable with thermal conductivity alone. Therefore, 
this is strong evidence that the Virgo cluster is being heated 
in the central regions by a mechanism other than thermal 
conduction. 

Despite the evidence for heating we are unable to make 
any strict quantitative predictions about the level of ther- 
mal conduction in Virgo. However, according to the observa- 
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Figure 6. Emissivity profiles for (from left to right) simulation 0k and 100k. The thick line shows the observed emissivity given by 
Ghizzardi et al. (2004), the thin solid line shows the emissivity after 6.34 X 10 8 yrs; the thin dashed line shows the emssivity at the end 
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Figure 7. Emissivity profiles for (from left to right) simulation 0.1k and Ik. The definitions for each line are as in figure l6l 



tional data 7 e g appears to vary significantly throughout the 
cluster. Our results suggest that large values of the thermal 
conductvity tend to smooth out variations of 7 e ff. In addi- 
tion, since thermal conductivity acts to reduce the tempera- 
ture gradient, any heating source located near the centre of 
a cluster in which thermal conduction operates at near the 
Spitzer value is likely to result in a much flatter temperature 
profile than is currently observed. In the case where thermal 
conduction is heavily suppressed, it is possible that an ad- 
ditional heat source could maintain the temperature in the 
cluster centre at the observed level while allowing the region 
outside to continue to cool radiatively and thus produce the 



observed wide bowl shape. Thus although we are unable to 
put strict constraints on the value of thermal conductivity 
in Virgo, our results suggest that it is probably s uppressed 
by a f actor of at least 10 or more. Recent work bv lBruggenl 
( 2003) and lRuszkowski fc Beeelmanl i2002l) suggests that, in 
1-d at least, it is possible to achieve a steady state between 
radiative cooling and simultaneous heating from an AGN 
and thermal conduction suppressed by factor similar to our 
suggestion. 
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Figure 9. Central density against central temperature evolving as a function of time for simulation 0.01k. The bottom right of the plot 
corresponds to the initial conditions and the system progresses upwards with passing time. 
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Figure 8. Comparison of the effective adiabatic index for simu- 
lation Ik (triple dot dash), 0k (dot dash) and the effective index 
from the fits to the observational data (solid line) and the data 
points. The triangles show the observational data that the fits are 
based on. At no point in time does the simulation data, within the 
central 15 kpc, approximate to the observational data. Therefore, 
for clarity, the 7 e g profiles for the simulation data are calculated 
for times where the central temperature is near its lowest point. 
The vertical line near the left of the figure is the error bar from 
the data point closest to the cluster centre. 



4.4 Central Temperature and density 

It is also instructive to plot the central density against cen- 
tral temperature, at various points in time throughout each 
simulation. Figure [5] shows this for simulation 0.01k. Each 
simulation shows the same three main characteristics: an ini- 



tial power-law dependence between falling temperature and 
increasing density followed by a turning point near which the 
temperature begins to increase up to an asymptotic value 
while the density continues to increase. The results for each 
individual simulation are given in table 2 where a and b are 
defined by a power-law relating temperature and density: 
n — bT a . Given that the early evolutionary behaviour of the 
central temperature and density is accurately described by 
a power-law it is also appropriate to describe this behaviour 
in terms of an effective adiabatic index, 7 e g which for the 
specified power-law is given by 1 + 1/a. The results given 
in table 2 show that the intial cooling phases are well de- 
scribed by power-laws and that a decreases for increasing 
values of thermal conductivity and b, which relates to the 
entropy, increases. 7 e g increases for larger thermal conduc- 
tivities because of the heating effect of thermal conduction 
in the cluster centre. 

After the initial cooling in which the temperature falls 
with time, we observe without exception a period in which 
the central temperature begins to rise to an asymptotic 
value of roughly 1.6 keV which appears to be the same for 
each simulation. We have been unable to investigate the be- 
haviour past this asymptotic point due to the smallness of 
the timestep at the time when this behaviour was observed. 

The temperature rise is probably due to the cooling 
catastrophe which occurs after 10 9 yrs. The rapid inflow 
of material in the cluster centre results in some compres- 
sional heating. After a certain time the central tempera- 
ture has reached its maximum and begins to decrease again. 
Throughout this behaviour the density continues to increase 
while the temperature probably oscillates around the asymp- 
totic value to some degree. 

It is possible that the central collapse causes some re- 
heating of the very centre of the cluster. Note that we did 
not observe this 'bounce' for the trial simulations conducted 
for l/8th of the volume, but with the same spatial resolution 
in the gas flow, of the simulations presented here. We also 
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find no evidence for a shock occuring near or after the turn- 
ing point. It appears that the power-law observed in these 
simulations and the other notable characteristics are unaf- 
fected by increasing the numerical resolution, see Appendix 

El 

4.5 Mass deposition rates 

In this Section we calculate mass flow rates as a function of 
radius. Thus we can directly estimate the rate at which gas 
flows in to the cluster centre for different thermal conduc- 
tivities. In addition, we also calculate theoretical mass flow 
rates and compare them to the results from our simulations. 

The mass deposition rate is given by the mass flux 
through a spherical surface 

= J pv.dA — 4nr 2 pv r [g s -1 ], (18) 

where v T is the projection of the velocity vector on to the 
radial vector. 

Figures 1101 to 1111 show the mass deposition rates as a 
function of radius for two different times. These figures show 
the local mass flow rates and not the integrated mass flow 
rate, M(<r). The first line (solid) is the mass deposition rate 
profile after 6.34 x 10 s yrs, the second line (dashed) is the 
profile for the final point in the simulation given in table 1. 

For the mass flow rate profiles after 6.34 x 10 s yrs the 
flow rate tends to zero in the cluster centre in all simulations. 
There is also a point near 150 kpc at which the flow rate 
reaches a maximu m. From here on we refer to this po int as 
the break-radius (lEdee et alJll992t lAllen et alJl200ll) . The 
increase outside the break-radius (at roughly 200 kpc) is 
probably not realistic, but is due to edge effects in which the 
initial velocities are large, but dissipate over the duration of 
the simulation. 

For the later profiles, when the cooling flow is well es- 
tablished, the central mass deposition rate increases up to 
several solar masses per year in all cases. We also note that 
the break-radius has moved further away from the cluster 
centre. The peak mass flow, at the break radius, is smaller 
for larger values of thermal conductivity. This is possibly 
a consequence of the action of shear viscosity which slows 
the inflow of the gas. In addition, the width of the mass 
flow rate peak increases with increasing thermal conductiv- 
ity and viscosity, indicating that in such cases there is more 
mass flowing towards the centre and that the inflow velocity 
falls off less quickly with radius. Thus it is possible that the 
central galaxies in clusters, which are thermally conducting, 
may be more massive, by a factor of a few, than those in 
cluster in which thermal conduction does not occur. 

The calculated mass deposition rate s are import ant for 
comparison with observational results. lEded (1200 ll) finds 
the mass of molcular gas around M87 to be less than 
1.3 x 1O 8 M0 . If the Virgo cluster has been cooling for roughly 
10 Gyr then the average mass deposition rate over this pe- 
riod cannot be greater than 0.013 Mq yr _1 . We find that for 
a thermal conductivity at the Spitzer value the mass deposi- 
tion rate in the centre of the cluster is negligible for the first 
~ 10 9 years, but rises to several Mq yr _1 thereafter. In ad- 
dition, the rate of mass deposition will continue to increase 
after this time as the cooling catastrophe accelerates. For 



the cases with suppressed thermal conductivity the central 
mass deposition rate is significant from earlier times mean- 
ing that more gas is deposited in the cluster centre for these 
simulations. Therefore, although the central mass flow rates 
of a few M©yr _1 appear at first rather modest, it is worth 
noting that the mass flow rates in the centre will quickly 
exceed the observations' upper limit. 

It is possible to compare the calculated mass flow rates 
from our simulations with analytical predictions. In order 
to do this we require both the density and velocity profiles. 
However, it is not possible to solve the full hydrodynamic 
equations analytically for the velocity profile as a function 
of time due to the complexity of the problem. Instead we 
use both observational and simulation data to constrain the 
functional form of the velocity profile at any given time. 

Figures 1101 and 1111 show that there is clearly a break 
radius in the m ass flow rat e profi les w hich agrees with th e 
observations of lEdee et all (Il992ft and lAllen et alJ fcOOll) . 
However, unlike these authors we plot the mass flow rate 
at each radius and not the cumulative mass flow rate. Since 
the radiative cooling rate is greatest near the cluster centre, 
gas in the central regions will lose its pressure support. This 
means that not only would we expect the gas infall velocity 
to tend asymptotically to zero for large radii, but that the 
scale height of the velocity profile will increase as a function 
of time as the less dense gas at larger radii has time to cool. 
The radiative cooling time as a function of radius is given by 
equation J7J . We assume the cluster gas to be isothermal and 
the density to follow a generalised /3-profile. Rearranging for 
the radius, we arrive at the cooling radius as a function of 
time 

'"cool 

/ / 2iA rado no \ w ,\ 2 ^ 
ro ~ U 3& bT°' 5 J J ' ( ' 

where ro is the scale height of the density distribution. For 
the Virgo cluster the value of /3 is roughly 0.47 for radii 
greater than roughly 30 kpc. 

Equation I19H implies that after times of order 0.3 Gyrs 
the velocity scale height grows proportional to f 1 /' 3 ' 3 '. 

Further constraints on the velocity profile can be deter- 
mined using equation 1181 and what we already know about 
the gas density of the Virgo cluster. The electron number 
density given by equation 110H implies that at large radii 
the density scales as r 1 ' 4 . Substituting this into equation 
1181 implies that in order for there to be a break radius 
the infall velocity of the gas must fall off more quickly than 
approximately r -0,6 . Furthermore, by the observation that 
the mass flow rate tends asymptotically to zero in the clus- 
ter centre (before the cooling catastrophe occurs) it follows 
that, at least, in the central region the velocity cannot fall off 
faster than r~ 2 . However, it should be noted that after the 
occurance of the cooling catastrophe the central mass flow 
rate becomes finite. This suggests that the velocity profile 
steepens at least in the central few kiloparsecs and may be 
an indication that the gas is in free-fall after the onset of the 
cooling catastrophe. Thus, for the Virgo cluster, if the infall 
velocity follows a power-law (v ~ r~ a ) then a may be con- 
strained to lie between 0.6 and 2. In general, for any cluster 
in which the gas density may be described by a /3-profile the 
constraints on a are (2-3/3)^Ca^C2. 

In the absence of an analytical expression for the veloc- 
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Table 2. Summary of the parameters describing the behaviour of the central temperature 
and density for each simulation. The temperatures are given in keV and densities in units of 
10 3 cm" 3 



k/k s 


a 


b 


7off 


turning point(T,n) 


asymptote(T) 





-1.60±0.025 


-30.7±0.41 


0.375±0.006 


~0.84, ~151 


~1.6 


0.01 


-1.67±0.016 


-29.7±0.27 


0.40±0.004 


~0.89, ~149 


~1.6 


0.1 


-1.96±0.119 


-24.5±1.95 


0.49±0.029 


~0.94, ~212 


~1.6 


1 


-2.32±0.070 


-18.4±1.17 


0.57±0.017 


~0.82, ~270 


~1.6 




Figure 10. Mass deposition profiles for (from left to right) simulations Ok and 0.01k. In all plots the solid line shows the mass flow rate 
at 6.34 X 10 s yrs; the dashed line refers to the end point in the simulation. The times of the end points are given in table 1. 
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ity profile we adopt a velocity profile as v ~ exp(— r 2 /r s 2 ) 
which satifies the condition above, where r s is the scale 
height of the velocity profile which is equivalent to the cool- 
ing radius described in equation 119H . 

In order to understand the dependence of the mass flow 
rate on velocity we use the density profile that we calculated 
for the initial conditions of these simulations and the veloc- 
ity profile described above and vary only the velocity scale 
height with time. The peak mass flow rate is determined by 
the velocity in the cluster centre, which by comparison with 
the simulation we chose to be of order 10 5 cm s -1 . The den- 
sity profile can be assumed to remain fixed since the den- 
sity should not change significantly until a cooling catas- 
trophe occurs. However, a gradual increase in the central 
density will inevitably occur due to the inflow of gas to- 
wards the cluster centre. The results of this basic analysis 
are compared to the simulation data in figure [HI It is ob- 
vious that despite the simplistic nature of the analysis its 
predictions match the simulation results reasonably well. In 
addition, the analysis should hold for clusters in general and 
not just the Virgo cluster. Thus our inflow velocity profile, 
v ~ exp(— r 2 /r s 2 ), can be taken as a generic approximation 
for the gas velocity in radiatively cooling clusters. 



5 SUMMARY 

Our results suggest that thermal conduction is unable to 
prevent the occurance of a cooling catastrophe and can only 
increase the cooling time in the cluster centre by a factor of 
a few. Tests using different coordinate systems and spatial 
resolutions show that the developement of a cooling catas- 
trophe is inevitable. We do not consider this to be a con- 
sequence of our chosen initial conditions since we allow the 
cluster to evolve self-consistently, taking into account all rel- 
evant physics of which we are aware, using standard intial 
conditions used in cluster formation theory. 

The temperature and density profiles from our simu- 
lated clusters fail to simultaneously converge to the observa- 
tional constraints. The mass flow rates we calculate also sug- 
gest that significantly more mass is deposited in the cluster 
centre than is observed. Somewhat contrary to expectations 
we also find that the mass flow rates are also greatest for 
the simulations in which thermal conductivion is most im- 
portant. The extra central concentration of mass this causes 
results in emissivities with broader peaks than for lower ther- 
mal conductivities. A consequence of this is that for these 
cases more energy has to be injected, for example, by a cen- 
tral AGN than for more weakly conducting cases to avert a 
cooling catastrophe. 

These results are in agreement with IVoigt et al.1 (120021) 
who found that in their sample thermal conduction could 
only balance radiative losses in the outer regions of clus- 
ter cor es. There is a slight c ontrad iction with the find- 
ings of Zakam ska fc Naravanl J2003ft who found that the 
temperature and density profiles of several galaxy clusters 
were consistent with a steady state maintained by phys- 
ically ineairhigjul_vahafis_rf therm al conduction. However, 
since IZakamska fc Naravanl ^2003) do not take into account 
the effect of line cooling, which is more significant at the 
lower temperatures near cluster centres, and only look for 
the temperature and density profiles which will ensure a 



steady state, it is unclear whether t his study and our s are 
comparable. Furthermore, ne ither of IVoigt et alJ J2002I) nor 
IZakamska fc Naravanl i2003T) have made predictions about 
the Virgo cluster. It would therefore be of interest to un- 
dertake a similar stu d y, to this, for a cluster for which 
IZakamska fc Naravanl J2003) suggest thermal conduction 
could maintain a steady state. 

The simulated and observed effective adiabatic index 
indicates that there is a heat source at the centre of the 
Virgo cluster. In addition, given this and the temperature 
and density profiles we suggest that thermal conductivity 
is suppressed to at most 1/100 to 1/10 of the full Spitzer 
value. However, even if thermal conductivity and viscosity 
are suppressed by a factor of 100, the thermal conduction 
and viscous dissipation of sound waves generated in the ICM 
could be significant. 

For the overall energy budget we suggest that thermal 
conduction is unlikely to be able to suppress the cooling in 
the majority of clusters. The reason for this is simple: the 
physics underlying thermal conduction and radiative cooling 
is too different for the two processes to balance each other. 
The fact that these energy terms cannot balance means that 
the temperature and density profiles are likely to be unstable 
on cosmologically relevant timescales. Thermal conductivity 
cannot be responsible for the universa l temp erature profiles 
of clusters as observed bv I Allen et alJ i200ll) . 
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Figure 12. Comparison of the mass flow rates for increasing time for simulation 0.01k data (left) and from simple theoretical predictions 
at the same times. 
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APPENDIX A: NUMERICAL RESOLUTION 
CONVERGENCE 

The numerical convergence of the simulations was tested by 
running simulation 0.01ft: with varied spatial resolution. We 
compare the evolution of temperature and density profiles 
for otherwise identical simulations in which the minimum 
refinement, maximum refinement and region of increased 
initial resolution in the central region are all increased by 
one level of refinement. For the size of the computational 
box we employ for these simulations, a refinement level of 9 
corresponds to a resolution of 0.31 kpc; a refinement level of 
10 corresponds to 0.15 kpc. The initial enhanced resolution 
of the central region (refinement of 6) corresponds to 2.53 
kpc in the central ~16 kpc. We also investigate the effect of 
increasing this resolution by one level of refinement to a res- 
olution of 1.26 kpc. The resolution in the outer regions of the 
cluster is determined by the minimum refinement level. This 
is initially set to 3 (corresponding to roughly 20 kpc) , but is 
also tested at level 4 (roughly 10 kpc). The refinement pa- 
rameters corresponding to each simulation are summarised 
in table Al. 

From figure |A"T1 it is evident that only for a central res- 
olution of refinement level 6, or lower, do we observe sig- 
nificantly different results to any other simulation. In this 
particular case it appears that the central resolution is not 
sufficient for the inner regions to be in hydrostatic equilib- 
rium when the simulation is initialised. This results in a 
slight collapse of material near the cluster centre and there- 
fore significant compressional heating. For the other simu- 
lations the differences in the temperature and density plots 
are minimal. The central effective adiabatic index is also un- 



16 E.C.D. Pope, G. Pavlovski, C.R. Kaiser, H. Fangohr 




r i i ■ , i I ■ i 



33 4Q 60 <U tm 




i ID i DO 



























| • s 






F IP 






i 










0L0 







m 4a 6o aa rCC 





a <n jn jj rod 



in i co 
d^lonce [Kpc] 




™i 1 1 i i m 




iu 1' i 1 1 id J 



1 ID 1DO 



Figure Al. Comparison of simulation 0.01k temperature and density profiles at each of the different numerical resolutions at intervals 
of 3.16 X 10 8 yrs starting at 3.16 X 10 8 yrs. The results for test 1 are shown by dot dashed lines, test 2: dashed, test 3: triple dot dash, 
test 4: long dash, test 5: dot 
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Table Al. Summary of the refinement parameters of the numerical resolution tests. 
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Test 


Central Refinement 


Maximum refinement 


Minimum refinement 


1 


7 


9 


3 


2 


6 


9 


3 


3 


7 


10 


3 


4 


8 


9 


3 


5 


7 


9 


4 



changed and in good agreement across the range of different 
resolutions. This series of tests demonstrates that for a cen- 
tral refinement of level 7, a maximum refinement of level 9 
and a minimum refinement of level 3 the results are in close 
agreement with those achieved with higher resolutions. 



APPENDIX B: HYDRODYNAMIC 
DIFFERENCES BETWEEN 3-D AND 1-D 
COORDINATE SYSTEMS GEOMETRIES 

In this section we outline the differences between a purely 
spherically symmetric geometry and a full 3-d geometry, for 
our problem, with the aid of the hydrodynamic equations 
and trial simulations of the Virgo cluster. For these pur- 
poses we have re-done two of the simulations from the main 
paper with i) zero thermal conductivity and viscosity and 
ii) full Spitzer thermal conductivity and viscosity, in 1-d 
spherical coordinates, both at the same spatial resolution as 
our 3-d runs and higher spatial resolutions. The minimum 
level of refinement for the 1-d simulation performed at the 
same spatial resolution as the 3-d runs was 7 correspond- 
ing to roughly 1.27 kpc. For the higher spatial resolution 
1-d test the minimum refinement was level 9, correspond- 
ing to a spatial scale of nearly 0.32 kpc. In both cases the 
maximum level of the refinement was 12 corresponding to 
approximately 0.04 kpc. 

The results of these 1-d spherically symmetric simula- 
tions explicitly show that a cooling catastrophe always oc- 
curs in 1-d as well as 3-d. Furthermore, the cooling catastro- 
phe takes less time to evolve in the 1-d spherically symmetric 
case than the 3-d case and increasing spatial resolution has 
almost no effect. For the simulation with zero thermal con- 
ductivity the cooling catastrophe occurs after approximately 
1 x 10 9 yrs-roughly 3 x 10 8 yrs before the occurance in 3-d. 
For the simulation with full Spitzer thermal conductivity the 
cooling catastrophe occurs after roughly 4.2 x 10 9 yrs which 
is approximately 0.5 x 10 9 yrs less than in the 3-d case (see 
figures [EH and inn ■ 

Figure lB"3l gives an indication of the non-radial flow in 
the central regions of our simulated clusters, in the absence 
of thermal conduction and viscosity. If the differences for 
simulations using different geometries are real physical ef- 
fects, rather than numerical errors, then we should also ex- 
pect to see differences in the equations of hydrodynamics 
for the 1-d and 3-d cases. Since viscosity enters into the mo- 
mentum equation it is an appropriate choice to investigate 
any differences which may arise due to dimensional effects. 



The equ ation for the v iscous force per unit volume is 
given by re.g. lSarazinlll98d) 

F vls = r)(v 2 v + ^VV.v\ (Bl) 

where rj is the bulk viscosity and v is fluid velocity. 

Using the standard vector identity we can substitute for 
the divergence term 

V 2 v = V(V.w) - V x (V x v) (B2) 

so that equation <BH becomes 

fvi S = 7 ?Qv 2 « + ivx(Vx«)j. (B3) 

Even without the evidence from simulations, from equa- 
tion <B3t alone it is evident that viscous processes are differ- 
ent in 3-d and 1-d cases due to shear effects. In 3-dimensions 
the term involving the vector products constitutes a mech- 
anism for the dissipation of momentum, whereas in 1-d all 
vector products are zero, by definition, and so this mecha- 
nism is suppressed. 

It seems that equation <B3t is able to account for some 
of the differences between the 1-d and 3-d simulations in 
which full Spitzer thermal conductivity and viscosity were 
present. However, there was also a significant difference in 
the time at which the cooling catastrophe occured in the 
simulations in which diffusion processes were not present. 
Studying the complete momentum equation allows us to in- 
vestigate these differences. 

^+(d.V)« = F--VP + i/( |v 2 i; + iv x (V xv)\ 
dt y ' p \3 3 J 

(B4) 

where F is the force per unit mass, P is the fluid pressure 
and v is the kinematic viscosity. 

Using a second standard vector calculus identity for the 
second term on the left hand side of equation 1B4II we find 

(v.V)v = iv(w.w) - v x (V x v), (B5) 

and substituting into equation <B4L the total momentum 
equation is 

^-+-V(v.v)-vxCVxv) = F--VP+v(-V 2 j)+-Vx(Vxti) 

(jt 2 p V 3 3 

(B6) 

Since vector products are zero in 1-dimension, equation 
<B6> reduces to 

dv l d(v 2 ) _ p IdP 4 1 d ( r 2dv\ , BJ) 
dt 2 dr p dr 3 r 2 dr \ dr J ' 
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Figure Bl. Temperature and density profiles evolving with time for simulation Oft for a 1-d spherically symmetric geometry. The profiles 
are plotted for the same times as the equivalent temperature and density profiles in section f4. II with the same spatial resolution, except 
for the final profile which is plotted at 1 X 10 9 yrs. 
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Figure B2. Temperature and density profiles evolving with time for simulation Ik for a 1-d spherically symmetric geometry. The profiles 
are plotted for the same times as the equivalent temperature and density profiles in section [4. II with the same spatial resolution, except 
for the final profile which is plotted at 4.2 X 10 9 yrs. 



in 1-d spherically symmetric coordinates. 

From equation iB7l it is clear that there is an additional 
dissipation term in the 3-d case which is not present in 1-d, 
even in the absence of viscosity. Thus, in agreement with 
our simulations, it is reasonable to expect the behaviour of 
a fluid to be different in 1-d compared to 3-d, even in the 
absence of viscosity. 

We note that the overall properties of the simulations, 
e.g. the temperature and density profiles, remain roughly 
the same in 1-d spherically symmetric and 3-d coordinates. 
Therefore, for simulations such as these it seems possible 
to gain a reasonable estimate of the main properties and 
timescales for the given circumstances at a lower compu- 
tational expense than for 3-d simulations. However, in or- 
der to gain a more complete understanding of the physical 



processes at work it is necessary to perform even these ba- 
sic simulations in 3-d due to the non-spherically symmetric 
process which may occur in both viscous and inviscid fluid 
dynamics. 
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Figure B3. Zoom on the centre of the simulation box. Velocity cut plane (velocity vectors shown as cones) highlights the total 
velocity vectors, demonstrating the essential non-radial features of the flow. We have included two sets of particle tracers (grey flux 
tubes) originating from the opposite corners (first set of 30 particle tracers has centre at (0.450,0.550,0.450) and the second set at 
(0.550,0.450,0.550)). The coordinates are shown in fraction of total box width (1. X 10 24 cm). 



